Réaliser des cartes avec R

Karine Pietropaoli, CMW CNRS

I Gestion et manipulation d’objets spatiaux avec sf et dplyr

II La cartographie thématique avec cartography

III La cartographie dynamique avec leaflet

Installation et chargement des extensions utiles

Création d’un nouveau projet Rstudio et d’un nouveau script
Installation des extensions nécessaires

install.packages("sf") # Gestion et manipulation d'objets spatiaux

install.packages("dplyr") # Manipulation de tableaux de données

install.packages("cartography") # Création de cartes thématiques

install.packages("leaflet") # Création de cartes web dynamiques

Chargement des extensions

library(sf)

library(dplyr) 

library(cartography)

library(leaflet)

I Gestion et manipulation d’objets spatiaux avec `sf` et `dplyr`

Import d’objets spatiaux
Gestion des projections cartographiques
Manipulation des données avec dyplr
Affichage des objets spatiaux
Transformation des objets spatiaux

1. Importer de données de type shapefile

Rappel :

Le format propriétaire shapefile (extension .shp) est parmi les formats de fichiers d’images vectorielles les plus utilisés. Ce fichier .shp stocke toute l’information liée à la géométrie des objets qu’il décrit et s’accompagne toujours d’au moins deux autres fichiers (.dbf et .shx) qui contiennent des informations complémentaires sur ces objets

Données mobilisées et sources :

La fonction st_read permet d’importer les shapefiles (et les autres formats SIG) dans R.

  com <- st_read(dsn = "SHP/COMMUNE_IGN", layer = "COMMUNE")
  dep <- st_read(dsn = "SHP/DEPARTEMENT_IGN", layer = "DEPARTEMENT")
  grandlyon <- st_read(dsn = "SHP/COMMUNE_GLYON", layer = "adr_voie_lieu.adrcommune")

Nature des objets crées

class(com)

## [1] "sf"         "data.frame"

com est à la fois un objet sf et data.frame. Les objets sfsont des data.frames enrichis d’une colonne contenant des géométries (listes de coordonnées définissant des points, des lignes ou des surfaces).

2. Gestion des projections cartographiques

La fonction st_crs permet de connaître la projection de nos objets spatiaux

st_crs(com) # Voir la projection

## Coordinate Reference System:
##   No EPSG code
##   proj4string: "+proj=lcc +lat_1=44 +lat_2=49 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs"

st_crs(grandlyon)

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

  st_crs(com) == st_crs(grandlyon)

## [1] FALSE

Ici, nos données spatiales sont dans des projections différentes

Harmonisation de la projection de nos différents objets avec la fonction st_transform :

com <- st_transform(com, crs = 4326)
dep <- st_transform(dep, crs = 4326)

st_crs(com) # Voir la projection

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

3. Manipulation des données avec `dyplr`

les objets sf étant aussi des data.frames “classiques”, nous pouvons les manipuler comme tels et c’est tout l’intérêt !

names(com)

##  [1] "ID_GEOFLA"  "CODE_COM"   "INSEE_COM"  "NOM_COM"   
##  [5] "STATUT"     "X_CHF_LIEU" "Y_CHF_LIEU" "X_CENTROID"
##  [9] "Y_CENTROID" "Z_MOYEN"    "SUPERFICIE" "POPULATION"
## [13] "CODE_ARR"   "CODE_DEPT"  "NOM_DEPT"   "CODE_REG"  
## [17] "NOM_REG"    "geometry"

Par exemple, pour travailler uniquement sur la région Auvergne Rhone-Alpes :

com_RA <- com %>% filter(CODE_REG == "84")
dep_RA <- dep %>% filter(CODE_REG == "84")

ou sur la ville de Lyon :

code <- 69381:69389
arr_lyon <- com_RA  %>% filter(INSEE_COM %in% code )

Pour ajouter des données à un objet sf, la jointure entre nos deux objets se fait comme pour deux data-frames classiques :

load("base_insee.Rdata") # chargement des données 

com_RA <- left_join(com_RA, base_insee, by = c("INSEE_COM"= "CODGEO"))

source des données : base comparateur de territoires INSEE

4. Affichage des objets

la fonction st-geometryassociée à la fonction de base plot permet un affichage simple des objets

st_geometry(dep_RA) %>% plot()

Superposition des objets

st_geometry(grandlyon) %>% plot(lwd = 0.5, col = NA, border = "grey40")
st_geometry(arr_lyon) %>% plot(lwd = 0.7, col= NA, border="red", add=TRUE)

5. Transformation des objets

Agrégation de polygones avec la fonction st_union

lyon  <- st_union(arr_lyon)

# résultat
st_geometry(lyon) %>% plot()

Extraire les centroïdes des polygones avec la fonction st_centroid

arr_lyon_c <-  st_centroid(arr_lyon)
st_geometry(arr_lyon) %>% plot()
st_geometry(arr_lyon_c) %>% plot(pch=19, col="red", add=T)

II Cartographie thématique avec `cartography`

Carte choroplèthe
Carte en symbole proportionnel

1. Carte choroplèthe

Rappel définition :

Une carte choroplèthe (du grec χώρος (khorê) : « zone/région » et πληθαίν (plethos) : « multiple ») est une représentation de quantités relatives à des espaces ou aires géographiques (tels la densité de population ou le revenu par habitant) par le moyen d’une « échelle » de tons gradués.

Calcul densité de population

com_RA$dens <- 100 * com_RA$POPULATION/com_RA$SUPERFICIE # superficie en hectare 
# 1 km2 = 100 hectares

# principaux indicateurs statistiques de la variable `dens`
summary(com_RA$dens)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##     0.079    17.324    48.346   179.948   120.896 19158.824

Carte choroplèthe avec la fonction chorolayer

choroLayer(com_RA, var = "dens")

Nombre de classes et discrétisation paramétrables

# nombre de classes et discrétisation paramétrables >>  9 méthodes disponibles
# The "q6" method uses the following quantile probabilities: 
#0, 0.05, 0.275, 0.5, 0.725, 0.95, 1.
choroLayer(com_RA, var = "dens" , method="q6" , nclass=6 )

“Habillage”

# Défintion d'une palette de couleurs avec la fonction `carto.pal`
mypalette <- carto.pal("orange.pal", 6) 

#   Contrôle de la légende et bordure 
choroLayer(com_RA, var = "dens" , col=mypalette,  border=NA, nclass = 6,
method="q6", legend.pos="topright", legend.title.txt = "")

# Ajout du nom des départements avec la fonction `labelLayer`
labelLayer(dep_RA, txt = "NOM_DEPT", cex = 0.7)

# Habillage titre, source, auteur, flêche nord... avec la fonction `layoutLayer`
layoutLayer(title = "Densité de population\n(habitants/km²)", sources = "GEOFLA® 2.2 IGN", 
  north = TRUE, author = "KP@2018" , scale = NULL,  col = NA , coltitle = "black")

# Ajout de la limite des départements
st_geometry(dep) %>% plot(lwd = 0.7, col = NA, border = "black", add = TRUE)
st_geometry(dep_RA) %>% plot(lwd = 0.7, col = NA, border = "black", add = TRUE )

2. Carte en symboles proportionnels

# extrait les communes du département de la Savoie
savoie <- com_RA %>% filter(DEP=="73")

# principaux indicateurs statistiques de la variable `P14_POP` (population 2014)
summary(savoie$P14_POP)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    41.0   338.2   578.0  1462.1  1161.0 59490.0

Carte en cercles proportionnels avec la fonction propSymbolsLayer

# affichage des communes du département
st_geometry(savoie) %>% plot(col = "#F1EEE8", border = "grey40")
# affichage des cercles proportionnels avec la fonction `propSymbolsLayer`
propSymbolsLayer(savoie, var = "P14_POP")

“Habillage”

# affichage des communes du département et contrôle des couleurs
st_geometry(savoie) %>% plot(col = "lightblue", border = "lightblue4", bg="lightblue1")

# contrôle des couleurs des symboles et de  la légende 
propSymbolsLayer(savoie, var = "P14_POP", 
                  col =  "white",
                  border = "grey",
                 legend.title.txt = "",
                 legend.style = "e")


# Habillage titre, source, auteur, flêche nord... avec la fonction `layoutLayer`
layoutLayer(title = "Nombre d'habitants dans le département de la Savoie",
            sources = "INSEE", theme = "blue.pal",
            scale = NULL, north = TRUE , author = "KP@2018")

# Ajout de la limite des départements alentour
st_geometry(dep_RA) %>% plot(lwd = 1, col = NA, border = "grey40", add = TRUE)

III Cartographie dynamique avec `leaflet`

Tout ce que vous avez tjrs voulu savoir sur `Leaflet`…

extension developpée par RStudio
permet l’utilisation de la bibliothèque Java-script libre de catographie en ligne du même nom
S’appuie sur les fonds de carte libre d’OpenStreetMap
S’appuie sur le WGS 84 (World Geodetic System 1984 : système géodésique mondial, révision de 1984) utilisé notamment par le système GPS.

Initialisation d’une carte et appel à OpenStreetMap

leaflet()   %>% # initialisation de la carte 
addTiles() #  ajout des carreaux OpenStreetMap par défaut

Ajout d’une couche `addCircleMarkers`

leaflet()   %>% # initialisation de la carte 
addTiles()  %>% #  ajout des carreaux OpenStreetMap par défaut
addCircleMarkers(data= arr_lyon_c)  # ajout d'un un point

Ajout d’une couche `addPolygons`

leaflet()   %>% # initialisation de la carte 
addTiles()  %>% #  ajout des carreaux OpenStreetMap par défaut
addCircleMarkers(data= arr_lyon_c) %>%  # ajout d'un point 
addPolygons(data=arr_lyon) # ajout d'un polygone

Des options graphiques

leaflet() %>% 
addTiles() %>%
addCircleMarkers(data= arr_lyon_c,# ajout des cercles
radius = 4,   # contrôle de la taille du cercle
fillColor = "brown", fillOpacity = 1,# gestion de la couleur et de l'opacité de l'intérieur du cercle
stroke = F) %>%  # élimine le trait du cercle
addPolygons(data=arr_lyon, # ajout polygone
popup = arr_lyon$NOM_COM,# affichage d'un label au clic
color = "brown", opacity = 0.5 , weight=2, fillColor = "brown")  # gestion de la couleur, de l'opacité et de l'épaisseur du trait

Carte en cercles proportionnels

leaflet() %>% 
addTiles() %>% 
addPolygons(data=arr_lyon) %>% 
addCircleMarkers(data= arr_lyon_c, 
fillColor = "brown", fillOpacity = 1, stroke = FALSE,
radius = sqrt((arr_lyon_c$POPULATION/100)/3.14), # aire cercle = pi * R^2 
popup= as.character(arr_lyon_c$POPULATION))

Cartes choroplèthes

Utilisation des palettes de couleurs de l’extension RColorBrewer

# installation et chargement de l'extension
install.packages ("RColorBrewer")

library(RColorBrewer)

# visualisation des différentes palettes de RColorBrewer
display.brewer.all()

Préparation des données

# filtre les communes de la Haute-Savoie
hsavoie <- com_RA %>% filter (DEP=="74")

# résumé statistique de la variable "MED14", (revenu médian 2014)
summary(hsavoie$MED14)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   18881   23114   25651   26676   29251   45902       1

variable continue : `colorNumeric`

# définition palette
palette <- colorNumeric( palette = "YlOrRd", domain = hsavoie$MED14)
leaflet() %>% 
addTiles() %>% 
addPolygons(data=hsavoie, color="white", weight=2, fillColor = palette(hsavoie$MED14),  
popup= as.character(round(hsavoie$MED14,0)), fillOpacity=1) %>% 
addLegend(pal = palette, values = hsavoie$MED14, position="topright", opacity=1,
title="Niveau de vie médian") # ajout de la légende

variable discrétisée :`colorBin`

cl <- c(18000, 23000, 28000, 33000, 38000, 43000, Inf) # définition des bornes de classes
# définition palette
palette <- colorBin( palette = "YlOrRd", domain = hsavoie$MED14, bins = cl) 
# carte
leaflet() %>% 
addTiles() %>% 
addPolygons(data=hsavoie, color="white", weight=2, fillColor = palette(hsavoie$MED14),  
popup= as.character(round(hsavoie$MED14, 0)), fillOpacity=1) %>% 
addLegend(pal = palette, values = hsavoie$MED14, position="topright", opacity=1,
title="Niveau de vie médian")

variable discrétisée : `colorQuantile`

# définition palette
palette <- colorQuantile( palette = "YlOrRd", domain = hsavoie$MED14) 
# carte
leaflet() %>% 
addTiles() %>% 
addPolygons(data=hsavoie, color="white", weight=2, fillColor = palette(hsavoie$MED14),  
popup= as.character(round(hsavoie$MED14, 0)), fillOpacity=0.6) %>% 
addLegend(pal = palette, values = hsavoie$MED14, position="topright", opacity=0.6,
title="Quartile niveau de vie")

Des options “interactives”

palette <- colorNumeric(palette = "YlOrRd", domain = hsavoie$MED14)
leaflet() %>% 
addTiles() %>% 
addPolygons(data=hsavoie, color="white", weight=2, fillColor = palette(hsavoie$MED14),  
label = paste(hsavoie$NOM_COM, "-", round(hsavoie$MED14,0)), fillOpacity=1, 
highlight = highlightOptions(bringToFront = TRUE, color = "grey")) %>% 
addLegend(pal = palette, values = hsavoie$MED14, position="topright", opacity=1,
title="Niveau de vie médian")